#include "common.cuh"

// A, B, C are column major
// X is along a column
// Y is along a row




__global__ void SGEMM(_float* A, _float* B, _float* C, int heightA, int widthA)
{
	const int heightB = widthA;
	const int heightC = heightA;
	const int bx = blockIdx.x; const int by = blockIdx.y;
	const int tx = threadIdx.x; const int ty = threadIdx.y;
	// compute pointers in A, B, C using thread ID
	_float* pA = A + VEC_LENGTH * bx + BLOCK_X * ty + tx;
	_float* pB = B + heightB *( BY * by + ty) + tx;
	_float* pC = C + heightC * (BY * by) + BLOCK_X * ty + VEC_LENGTH * bx + tx;

	_float* pBEnd = pB + widthA;
	__shared__ _float sharedB[BX][BY+1];
	_float c[BY] = {0, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0 ,0 ,0 , 0 ,0 ,0, 0};
	
	while(pB < pBEnd)
	{
//#pragma unroll
		for(int i=0; i < BY; i+=4)
			sharedB[tx][ty+i] = pB[i*heightB];
		pB += BX;
		__syncthreads();
//#pragma unroll
		for(int i=0; i < BY; ++i)
		{
			c[0] += pA[0] * sharedB[i][0];
			c[1] += pA[0] * sharedB[i][1];
			c[2] += pA[0] * sharedB[i][2];
			c[3] += pA[0] * sharedB[i][3];
			c[4] += pA[0] * sharedB[i][4];
			c[5] += pA[0] * sharedB[i][5];
			c[6] += pA[0] * sharedB[i][6];
			c[7] += pA[0] * sharedB[i][7];
			c[8] += pA[0] * sharedB[i][8];
			c[9] += pA[0] * sharedB[i][9];
			c[10] += pA[0] * sharedB[i][10];
			c[11] += pA[0] * sharedB[i][11];
			c[12] += pA[0] * sharedB[i][12];
			c[13] += pA[0] * sharedB[i][13];
			c[14] += pA[0] * sharedB[i][14];
			c[15] += pA[0] * sharedB[i][15];

			pA += heightA;
		}
		__syncthreads();
	}
	for(int i=0; i < BY; ++i, pC += heightC)
		pC[0] = c[i];
}